function [gx,hx,exitflag,gxhx_info]=gx_hx(fir_ord,stake)

%[gx,hx,exitflag]=gx_hx(fy,fx,fyp,fxp,stake);
%computes the matrices gx and hx that define the first-order approximation to the solution  
%of a dynamic stochastic general equilibrium model. 
%Following the notation in Schmitt-Grohe and Uribe (JEDC, 2004), the model's equilibrium conditions 
%take the form
%E_t[f(yp,y,xp,x)=0.
%The solution is of the form
%xp = h(x,sigma) + sigma * eta * ep
%y = g(x,sigma).
%The first-order approximations to the functions g and h around the point (x,sigma)=(xbar,0), where xbar=h(xbar,0), are:
%h(x,sigma) = xbar + hx (x-xbar) 
%and
%g(x,sigma) = ybar + gx * (x-xbar),
%where ybar=g(xbar,0). 
%The variable exitflag takes the values 0 (no solution), 1 (unique solution), 2 (indeterminacy), or 3 (z11 is not invertible).
%Inputs: fy fyp fx fxp stake
% The parameter stake ensures that all eigenvalues of hx are less than stake in modulus (the default is stake=1).
%Outputs: gx hx exitflag
%This program is a modified version of solab.m by Paul Klein  (JEDC, 2000).
%Modified on August 30 2009 by Ryan Chahrour (rc2374@columbia.edu) to replace qzdiv.m by ordqz.m  to increase speed.
%Original Version by Stephanie Schmitt-Grohe and Martin Uribe
%Revised Version by Xing Guo
%Main Revision: make the computation more robust to potential rank
%               insufficiency
%Date July 17, 2001, May 11 2006, November 19, 2009

fy=fir_ord.YY;
fx=fir_ord.XX;
fyp=fir_ord.YYp;
fxp=fir_ord.XXp;
% stake=1;
if nargin<2
    stake=1;
end
exitflag = 1;

%Create system matrices
A = full([-fxp -fyp]);
B = full([fx fy]);
NK = size(fx,2);

%Complex Schur Decomposition
[s,t,q,z] = qz(A,B);   


%Pick non-explosive (stable) eigenvalues
slt = (abs(diag(t))<stake*abs(diag(s)));  
gxhx_info.eigval    =   sort(abs(diag(t))./abs(diag(s)));
nk=sum(slt);
gxhx_info.nk        =   nk;
gxhx_info.NK        =   NK;
%Reorder the system with stable eigs in upper-left
[s,t,q,z] = ordqz(s,t,q,z,slt);   

%Split up the results appropriately
z21 = z(nk+1:end,1:nk);
z11 = z(1:nk,1:nk);

s11 = s(1:nk,1:nk);
t11 = t(1:nk,1:nk);

gx      =   [];
hx      =   [];

if nk>NK
   % warning('The Equilibrium is Locally Indeterminate')
    exitflag=2;
elseif nk<NK
   % warning('No Local Equilibrium Exists')
    exitflag = 0;
else
    [V,D]   =   eig(z11);
    d       =   abs(diag(D));  
    tol     =   eps*1e4;
    if min(d)<tol && nk==NK
       % warning('Invertibility condition violated')
        exitflag    =   3;
        ind         =   find(d<tol);
        D(ind,ind)  =   D(ind,ind)+tol;
    end

    %Compute the Solution
    z11i    =   V/D/V;
    gx      =   real(z21*z11i);  
    hx      =   real(z11*(s11\t11)*z11i);
end


% switch exitflag
%     case 1 
%         z11i = z11\eye(nk);
%         gx = real(z21*z11i);  
%         hx = real(z11*(s11\t11)*z11i);
%     case 3
%         z22     =   z(nk+1:end,nk+1:end);
%         if rank(z22)==size(fy,2)
%             z12     =   z(1:nk,nk+1:end);
%             gx      =   real(-(z22')\z12');
%             M       =   z11'-z21'*gx;
%             if rank(M)==nk
%                 exitflag    =   1;
%                 hx  =   real(M\(s11\t11)*M);
%             end
%         end
% end
